Boring setup stuff

reqpkg <- c("magrittr","dplyr","ggplot2","ggpubr","DT","reshape2")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "magrittr"
## [1] '1.5'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
## [1] "DT"
## [1] '0.13'
## [1] "reshape2"
## [1] '1.4.3'
reqpkg <- c("purrr","nplr","car","compute.es","effects","multcomp","pastecs","WRS2","outliers")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "purrr"
## [1] '0.3.3'
## [1] "nplr"
## [1] '0.1.7'
## [1] "car"
## [1] '3.0.7'
## [1] "compute.es"
## [1] '0.2.4'
## [1] "effects"
## [1] '4.1.4'
## [1] "multcomp"
## [1] '1.4.12'
## [1] "pastecs"
## [1] '1.3.21'
## [1] "WRS2"
## [1] '1.0.0'
## [1] "outliers"
## [1] '0.14'
theme_set(theme_pubr())
select <- dplyr::select
colors <- c("#FC2A1C", "#103EFB")
uc_colors <- c("#800000", "#D6D6CE")

path <- "../data/R/Obesity/"
# Import data ####

cohort1_1 <- read.csv2(paste0(path, "1_pt1_total.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort1.1_",colnames(.)))
cohort1_1_dark <- read.csv2(paste0(path, "1_pt1_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort1.1_",colnames(.)))
cohort1_1_light <- read.csv2(paste0(path, "1_pt1_light.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort1.1_",colnames(.)))
cohort1_1_covar <- read.csv2(paste0(path, "1_pt1_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

# cohort1_2 <- read.csv2(paste0(path, "1_pt2_total.csv"), sep = ",") %>% 
#   set_colnames(paste0("cohort1.2_",colnames(.)))
# cohort1_2_dark <- read.csv2(paste0(path, "1_pt2_dark.csv"), sep = ",") %>% 
#   set_colnames(paste0("cohort1.2_",colnames(.)))
# cohort1_2_light <- read.csv2(paste0(path, "1_pt2_light.csv"), sep = ",") %>% 
#   set_colnames(paste0("cohort1.2_",colnames(.)))
# cohort1_2_covar <- read.csv2(paste0(path, "1_pt2_covar.csv"), sep = ",") %>% 
#   mutate_all(function(x) as.numeric(as.character(x)))

cohort1_3 <- read.csv2(paste0(path, "1_pt3_total.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort1.3_",colnames(.)))
cohort1_3_dark <- read.csv2(paste0(path, "1_pt3_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort1.3_",colnames(.)))
cohort1_3_light <- read.csv2(paste0(path, "1_pt3_light.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort1.3_",colnames(.)))
cohort1_3_covar <- read.csv2(paste0(path, "1_pt3_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

cohort2_1 <- read.csv2(paste0(path, "2_pt1_total.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.1_",colnames(.)))
cohort2_1_dark <- read.csv2(paste0(path, "2_pt1_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.1_",colnames(.)))
cohort2_1_light <- read.csv2(paste0(path, "2_pt1_light.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.1_",colnames(.)))
cohort2_1_covar <- read.csv2(paste0(path, "2_pt1_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

cohort2_2 <- read.csv2(paste0(path, "2_pt2_total.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.2_",colnames(.)))
cohort2_2_dark <- read.csv2(paste0(path, "2_pt2_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.2_",colnames(.)))
cohort2_2_light <- read.csv2(paste0(path, "2_pt2_light.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.2_",colnames(.)))
cohort2_2_covar <- read.csv2(paste0(path, "2_pt2_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

cohort2_3 <- read.csv2(paste0(path, "2_pt3_total.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.3_",colnames(.)))
cohort2_3_dark <- read.csv2(paste0(path, "2_pt3_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.3_",colnames(.)))
cohort2_3_light <- read.csv2(paste0(path, "2_pt3_light.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.3_",colnames(.)))
cohort2_3_covar <- read.csv2(paste0(path, "2_pt3_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

cohort3_1 <- read.csv2(paste0(path, "3_pt1_total.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.1_",colnames(.)))
cohort3_1_dark <- read.csv2(paste0(path, "3_pt1_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.1_",colnames(.)))
cohort3_1_light <- read.csv2(paste0(path, "3_pt1_light.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.1_",colnames(.)))
cohort3_1_covar <- read.csv2(paste0(path, "3_pt1_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

cohort3_2 <- read.csv2(paste0(path, "3_pt2_total.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.2_",colnames(.)))
cohort3_2_dark <- read.csv2(paste0(path, "3_pt2_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.2_",colnames(.)))
cohort3_2_light <- read.csv2(paste0(path, "3_pt2_light.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.2_",colnames(.)))
cohort3_2_covar <- read.csv2(paste0(path, "3_pt2_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

cohort3_3 <- read.csv2(paste0(path, "3_pt3_total.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.3_",colnames(.)))
cohort3_3_dark <- read.csv2(paste0(path, "3_pt3_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.3_",colnames(.)))
cohort3_3_light <- read.csv2(paste0(path, "3_pt3_light.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.3_",colnames(.)))
cohort3_3_covar <- read.csv2(paste0(path, "3_pt3_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

cohort2_run1 <- read.csv2(paste0(path, "2_macro13.csv_macro13.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort2.1_",colnames(.)))
cohort3_run1 <- read.csv2(paste0(path, "3_macro13.csv_macro13.csv"), sep = ",") %>% 
  set_colnames(paste0("cohort3.1_",colnames(.)))


covar_SPF_pt1 <- c(cohort1_1_covar[[1]][c(2,3)], cohort2_1_covar[[1]][c(6,7,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
covar_GF_pt1 <- c(cohort1_1_covar[[1]][1], cohort2_1_covar[[1]][c(1,2,3,4,5)], cohort3_1_covar[[1]][c(1,2,3)])
covar_SPF_pt3 <- c(cohort1_3_covar[[1]][c(2,3)], cohort2_3_covar[[1]][c(4,5,6,7)], cohort3_3_covar[[1]][c(4,5,6)])
covar_GF_pt3 <- c(cohort1_3_covar[[1]][1], cohort2_3_covar[[1]][c(1,2,3)], cohort3_3_covar[[1]][c(1,2,3)])

total2_pt1 <- cohort1_1 %>% 
  bind_cols(cohort2_1) %>% 
  bind_cols(cohort3_1) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

light_pt1 <- cohort1_1_light %>% 
  bind_cols(cohort2_1_light) %>% 
  bind_cols(cohort3_1_light) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

dark_pt1 <- cohort1_1_dark %>% 
  bind_cols(cohort2_1_dark) %>% 
  bind_cols(cohort3_1_dark) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

total_pt1 <- bind_rows(light_pt1, dark_pt1) %>% 
  mutate_at(dplyr::vars(contains("VO2")), ~ ifelse(. < 0.6, NA, .))

total2_pt2 <- cohort2_2 %>%
  bind_cols(cohort3_2) %>%
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

light_pt2 <- cohort2_2_light %>%
  bind_cols(cohort3_2_light) %>%
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

dark_pt2 <- cohort2_2_dark %>%
  bind_cols(cohort3_2_light) %>%
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

total_pt2 <- bind_rows(light_pt2, dark_pt2) %>% 
  mutate_at(dplyr::vars(contains("VO2")), ~ ifelse(. < 0.6, NA, .))

total2_pt3 <- cohort1_3 %>% 
  bind_cols(cohort2_3) %>% 
  bind_cols(cohort3_3) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

light_pt3 <- cohort1_3_light %>% 
  bind_cols(cohort2_3_light) %>% 
  bind_cols(cohort3_3_light) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

dark_pt3 <- cohort1_3_dark %>% 
  bind_cols(cohort2_3_dark) %>% 
  bind_cols(cohort3_3_dark) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

total_pt3 <- bind_rows(light_pt3, dark_pt3)

total_run <- cohort2_run1 %>% 
  bind_cols(cohort3_run1) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))
  

# Functions ####

select_GF <- function(df) {
  d1 <- df %>% select(contains(c("cohort1.1", "cohort1.2")))
  for (name in colnames(d1)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(6,9)) {
      d1 %<>% select(-all_of(name))
    }
  }
  
  d2 <- df %>% select(contains("cohort1.3"))
  for (name in colnames(d2)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(11,12)) {
      d2 %<>% select(-all_of(name))
    }
  }
  
  d3 <- df %>% select(contains(c("cohort2.1","cohort2.2")))
  for (name in colnames(d3)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(6,9,10,11)) {
      d3 %<>% select(-all_of(name))
    }
  }
  
  d4 <- df %>% select(contains("cohort2.3"))
  for (name in colnames(d4)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(9,10,11,12)) {
      d4 %<>% select(-all_of(name))
    }
  }
  
  d5 <- df %>% select(contains(c("cohort3.1","cohort3.2")))
  for (name in colnames(d5)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(12,13,14)) {
      d5 %<>% select(-all_of(name))
    }
  }
  
  d6 <- df %>% select(contains("cohort3.3"))
  for (name in colnames(d6)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(4,5,6)) {
      d6 %<>% select(-all_of(name))
    }
  }
  
  d <- bind_cols(d1, d2, d3, d4, d5, d6)
  return(d)
}

select_SPF <- function(df) {
  d1 <- df %>% select(contains(c("cohort1.1","cohort1.2")))
  for (name in colnames(d1)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(4,5)) {
      d1 %<>% select(-all_of(name))
    }
  }
  
  d2 <- df %>% select(contains("cohort1.3"))
  for (name in colnames(d2)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(10)) {
      d2 %<>% select(-all_of(name))
    }
  }
  
  d3 <- df %>% select(contains(c("cohort2.1","cohort2.2")))
  for (name in colnames(d3)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,2,3,4,5)) {
      d3 %<>% select(-all_of(name))
    }
  }
  
  d4 <- df %>% select(contains("cohort2.3"))
  for (name in colnames(d4)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,2,3)) {
      d4 %<>% select(-all_of(name))
    }
  }
  
  d5 <- df %>% select(contains(c("cohort3.1","cohort3.2")))
  for (name in colnames(d5)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(9,10,11)) {
      d5 %<>% select(-all_of(name))
    }
  }
  
  d6 <- df %>% select(contains("cohort3.3"))
  for (name in colnames(d6)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,2,3)) {
      d6 %<>% select(-all_of(name))
    }
  }
  
  d <- bind_cols(d1, d2, d3, d4, d5, d6)
  return(d)
}

graph_groups <- function(data, channel, title="", type="m", colors=c("#103EFB", "#FC2A1C")) {
  df <- data %>% select(contains(channel))
  df.observation <- 1:nrow(df)
  SPF <- df %>% select_SPF() %>% set_colnames(paste0("SPF_", colnames(.))) %>% melt(id.vars=NULL) %>% select(-variable)
  SPF$variable <- rep("SPF", nrow(SPF))
  GF <- df %>% select_GF() %>% set_colnames(paste0("GF_", colnames(.))) %>% melt(id.vars=NULL) %>% select(-variable)
  GF$variable <- rep("GF", nrow(GF))
  df <- bind_rows(SPF, GF)
  df$observation <- rep(df.observation, nrow(df)/length(df.observation))
  df <- na.omit(df)
  
  if (type=="m") {
    p <- df %>% group_by(variable, observation) %>% 
      summarise(value=mean(value)) %>% 
      ggplot(aes(x=observation, y=value, color=variable)) +
      geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
      geom_rect(aes(xmin=721, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
      geom_line(size=1) +
      scale_color_manual(values = colors) +
      scale_x_continuous(breaks=seq(0,nrow(df),240)) +
      ggtitle(title) +
      ylab(channel)
  } else if (type=="s") {
    p <- ggplot(df, aes(x=observation, y=value, color=variable)) + 
      geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
      geom_rect(aes(xmin=721, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
      scale_color_manual(values = colors) +
      scale_x_continuous(breaks=seq(0,nrow(df),240)) +
      stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", alpha = 0.7) +
      ggtitle(title) + 
      ylab(channel)
  }
  return(p)
}

graph_mice <- function(data, channel, title="") {
  df <- data %>% select(contains(channel))
  df.observation <- 1:nrow(df)
  SPF <- df %>% select_SPF() %>% set_colnames(paste0("SPF_", colnames(.))) %>% melt()
  GF <- df %>% select_GF() %>% set_colnames(paste0("GF_", colnames(.))) %>% melt()
  df <- bind_rows(SPF, GF) %>% set_colnames(c("mouse", "channel"))
  df$observation <- rep(df.observation, nrow(df)/length(df.observation))
  p <- ggplot(df, aes(x=observation, y=channel, color=mouse)) + 
    geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
    geom_rect(aes(xmin=721, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
    scale_x_continuous(breaks=seq(0,nrow(df),240)) +
    geom_line() +
    ggtitle(title) +
    ylab(channel)
  return(p)
}

signif.floor <- function(x, n){
  pow <- floor( log10( abs(x) ) ) + 1 - n
  y <- floor(x / 10 ^ pow) * 10^pow
  # handle the x = 0 case
  y[x==0] <- 0
  y
}

signif.ceiling <- function(x, n){
  pow <- floor( log10( abs(x) ) ) + 1 - n
  y <- ceiling(x / 10 ^ pow) * 10^pow
  # handle the x = 0 case
  y[x==0] <- 0
  y
}

find_range <- function(data, channel, segments=10) {
  data <- data %>% select(contains(channel)) %>% melt(id.vars = NULL)
  r <- range(data$value, na.rm=TRUE)
  cut <- max(data$value, na.rm=TRUE)-min(data$value, na.rm=TRUE)
  out <- list("range"=r, "cut"=cut/segments)
  return(out)
}


rcf <- function(data, channel, id, min, max, b=0.1) {
  #' Relative cumulative frequency function
  #'
  #' Takes a vector and returns a dataframe with relative cumulative frequency per break unit
  #'
  #' @param data numeric vector
  #' @param channel character. Channel to label the rcf output.
  #' @param id character. ID for this sample.
  #' @param min double. Minimum value of data.
  #' @param max double. Maximum value of data.
  #' @param b double. Number to cut vector.
  #' Before running, check for valid break values to cut the data. One way is to divide range of data by 10 `(max(data)-min(data))/10`.
  
  breaks <- seq(min, max, by=b)
  # breaks <- seq(signif.floor(min(data), 2), signif.ceiling(max(data), 2), by=b)
  duration.cut <- cut(data, breaks, right = F)
  duration.freq <- table(duration.cut)
  duration.cumfreq = cumsum(duration.freq) %>% prepend(0)
  duration.cumrelfreq = duration.cumfreq/length(data)
  out <- duration.cumrelfreq %>% as.data.frame()
  out <- cbind(breaks, out) %>% set_colnames(c("breaks", channel)) %>% set_rownames(NULL)
  out$id <- rep(id, nrow(out))
  return(out)
}

EC50 <- function(x) {
  out <- lapply(x, function(y) {
    nplr(y$breaks, y[, 2], silent = T, useLog = FALSE)
  }) %>% 
    sapply(function(z) {
      getEstimates(z, 0.5)$x
    })
  out %<>% set_rownames(NULL)
  return(out)
}

lm_eqn <- function(df, y, x){
  m <- lm(y ~ x, df);
  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
                   list(a = format(unname(coef(m)[1]), digits = 2),
                        b = format(unname(coef(m)[2]), digits = 2),
                        r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq));
}

lm_eqn_str <- function(df, y, x){
  m <- lm(y ~ x, df);
  eq <- sprintf("y = %s + %sx, r^2 = %s", format(unname(coef(m)[1]), digits = 2),format(unname(coef(m)[2]), digits = 2),format(summary(m)$r.squared, digits = 3))
  return(eq)
}

ANCOVA <- function(df) {
  # Takes dataframe with column 1 as independent variable (e.g., WT), column 2 as dependent variable, and column 3 as covariate
  
  cat("linear fit of dependent variable with covariates\n")
  p1 <- ggplot(df, aes_string(x=names(df)[2], y=names(df)[3])) +
    geom_smooth(method = "lm", se = FALSE) +
    geom_point() +
    # geom_text(x = mean(na.omit(df[[2]])), y = mean(na.omit(df[[3]])), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
    geom_text(x = mean(na.omit(df[[2]])), y = (mean(na.omit(df[[3]])) + IQR(df[[3]], na.rm = T)), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
    ggtitle("Linear model - DV and Covar") +
    xlab("DV") + ylab("Covar")
  print(p1)
  
  m <- lm(df[[3]] ~ df[[2]], df)
  
  summary(m) %>% print()
  lm_eqn_str(df, df[[3]], df[[2]]) %>% cat()
  
  plot1 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[2])) +
    geom_boxplot() +
    # ggtitle("DV vs IV") +
    xlab("IV") + ylab("DV")
  plot2 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[3])) +
    geom_boxplot() +
    # ggtitle("Covar vs IV") +
    xlab("IV") + ylab("Covar")
  p2 <- ggarrange(plot1, plot2,
                  labels = c("DV vs IV", "Covar vs IV"),
                  ncol = 2)
  # grid.arrange(plot1, plot2, ncol=2)
  print(p2)
  
  cat("\n―――――――――――――――――――――\n")
  cat("Some stats on the variables\n")
  cat("DV - IV:\n")
  by(df[[2]], df[[1]], stat.desc) %>% print()
  cat("Covar - IV:\n")
  by(df[[3]], df[[1]], stat.desc) %>% print()
  cat("levene's test\n")
  leveneTest(df[[2]], df[[1]], center = median) %>% print()
  
  cat("―――――――――――――――――――――\n")
  cat("ANOVA of DV with IV\n")
  model <- aov(df[[2]] ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("―――――――――――――――――――――\n")
  cat("ANOVA of covariate with IV, independent of DV\n")  
  model <- aov(df[[3]] ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("―――――――――――――――――――――\n")
  cat("ANOVA of DV with IV, normalized to covariate\n")
  c <- df[[2]]/df[[3]]
  model <- aov(c ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("―――――――――――――――――――――\n")
  cat("ANCOVA\n")
  model <- aov(df[[2]] ~ df[[1]] + df[[3]], data = df)
  summary(model) %>% print()
  out <- Anova(model, type= "III")
  print(out)
  return(out)
}

check_ANCOVA <- function(res) {
  if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] < 0.05) {
    print(paste("Varies significantly with DV and Covar,"), res$`Pr(>F)`[2], res$`Pr(>F)`[3], sep=" ")
  } else if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] > 0.05) {
    print(paste("Varies significantly with DV,", "p =", res$`Pr(>F)`[2], sep =" "))
  } else if (res$`Pr(>F)`[2] > 0.05 && res$`Pr(>F)`[3] < 0.05) {
    print(paste("Varies significantly with Covar,", "p =", res$`Pr(>F)`[3], sep=" "))
  }
}


rcf_pipe <- function(data, cycle, channel, segments = 50, id1 = "SPF", id2 = "GF", colors = c("#128D15", "#FD9226")) {
  r <- find_range(data, channel, segments)
  print("Range:")
  print(r)
  
  data <- data %>% select(contains(channel))
  d <- list()
  d[[id1]] <- data %>% select_SPF()
  d[[id2]] <- data %>% select_GF()
  
  rcf_res <- list()
  rcf_res[[id1]] <- lapply(d[[id1]], function(x) rcf(x, channel, id1, r$range[1], r$range[2], r$cut))
  rcf_res[[id2]] <- lapply(d[[id2]], function(x) rcf(x, channel, id2, r$range[1], r$range[2], r$cut))
  
  data <- bind_rows(bind_rows(rcf_res[[id1]]), bind_rows(rcf_res[[id2]])) %>% split(.$id)
  
  Models <- lapply(data, function(x){
    nplr(x$breaks, x[[channel]], silent = TRUE, useLog = F)
  })
  
  overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main=paste0("RCF of ", id1, " vs ", id2, " mice ", channel, ", ", cycle), cex.main=1.5, lwd = 3, Cols = c("#103EFB", "#FC2A1C"))
  abline(h=0.5)
  
  print(paste0(id1, " EC50: ", getEstimates(Models[[id1]], 0.5)$x))
  print(paste0(id2, " EC50: ", getEstimates(Models[[id2]], 0.5)$x))
  
  return(list(rcf_res[[id1]], rcf_res[[id2]]))
}

Summary table

Stage Cycle Variable SPF EC50 GF EC50 ANCOVA, DV ANCOVA, BM ANOVA, DV/BM ANOVA T test
1 Total VO2 1.750016179 1.608772088 0.008334 0.057229 0.000356 0.0568 0.05684
1 Light VO2 1.589060263 1.445741617 0.008267 0.076576 0.000387 0.0428 0.04277
1 Dark VO2 1.934114431 1.834036361 0.01839 0.06017 0.00113 0.123 0.1227
2 Total VO2 1.813006427 1.696862709 0.07927 0.17518 0.0197 0.219 0.2555
2 Light VO2 1.713093719 1.592411293 0.04997 0.16483 0.00653 0.142 0.1422
2 Dark VO2 1.921942202 1.840620907 0.191 0.2992 0.0785 0.367 0.3665
3 Total VO2 1.86394965 1.680559531 0.0002903 0.0051489 0.0163 0.704 0.01629
3 Light VO2 1.690084297 1.453820561 5.72E-05 0.003407 0.868 0.00352 0.003519
3 Dark VO2 2.05497862 1.933216996 0.005706 0.040283 0.504 0.0531 0.05315
1 Total RQ 0.8688879428 0.8263887997 0.03228 0.75872 0.000801 0.00477 0.00477
1 Light RQ 0.8083040868 0.7568620384 0.03168 0.78371 0.000622 0.00488 0.004882
1 Dark RQ 0.905887016 0.8881087771 0.1386 0.7982 0.0023 0.0888 0.08875
2 Total RQ 0.8899970429 0.8133868495 0.0106284 0.8573954 0.000897 0.00181 0.001961
2 Light RQ 0.8669847586 0.7878691449 0.0253839 0.5517931 0.00175 0.00346 0.003457
2 Dark RQ 0.9021469071 0.832689248 0.009808 0.849794 0.000779 0.00269 0.002689
3 Total RQ 0.8438159684 0.814793611 0.1087 0.5136 0.0267 0.0888 0.08881
3 Light RQ 0.7967254917 0.7508655951 0.1724047 0.6535429 0.0866 0.12 0.1196
3 Dark RQ 0.8717019528 0.866975323 0.545 0.6377 0.00795 0.689 0.6895
stage <- c("before HF diet", "immediately after HF diet", "1 month after HF diet")

SPF_total_VO2 <- c(1.750016179,1.813006427,1.86394965)
SPF_light_VO2 <- c(1.589060263,1.713093719,1.690084297)
SPF_dark_VO2 <- c(1.934114431,1.921942202,2.05497862)
GF_total_VO2 <- c(1.608772088,1.696862709,1.680559531)
GF_light_VO2 <- c(1.445741617,1.592411293,1.453820561)
GF_dark_VO2 <- c(1.834036361,1.840620907,1.933216996)

SPF_total_RQ <- c(0.8688879428,0.8899970429,0.8438159684)
SPF_light_RQ <- c(0.8083040868,0.8669847586,0.7967254917)
SPF_dark_RQ <- c(0.905887016,0.9021469071,0.8717019528)
GF_total_RQ <- c(0.8263887997,0.8133868495,0.814793611)
GF_light_RQ <- c(0.7568620384,0.7878691449,0.7508655951)
GF_dark_RQ <- c(0.8881087771,0.832689248,0.866975323)

c1 <- "VO2"
c2 <- "Total"

p1 <- data.frame(stage, SPF_total_VO2, GF_total_VO2) %>% melt() %>% 
  ggplot(aes(x=stage, y=value, group=variable, color=variable)) +
  geom_line() +
  ggtitle(paste(c2, "cycle", c1, "EC50", sep = " ")) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x=element_blank()) +
  scale_color_manual(values = colors)

c2 <- "Light"

p2 <- data.frame(stage, SPF_light_VO2, GF_light_VO2) %>% melt() %>% 
  ggplot(aes(x=stage, y=value, group=variable, color=variable)) +
  geom_line() +
  ggtitle(paste(c2, "cycle", c1, "EC50", sep = " ")) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x=element_blank()) +
  scale_color_manual(values = colors)

c2 <- "Dark"

p3 <- data.frame(stage, SPF_dark_VO2, GF_dark_VO2) %>% melt() %>% 
  ggplot(aes(x=stage, y=value, group=variable, color=variable)) +
  geom_line() +
  ggtitle(paste(c2, "cycle", c1, "EC50", sep = " ")) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x=element_blank()) +
  scale_color_manual(values = colors)

c1 <- "RQ"
c2 <- "Total"

p4 <- data.frame(stage, SPF_total_RQ, GF_total_RQ) %>% melt() %>% 
  ggplot(aes(x=stage, y=value, group=variable, color=variable)) +
  geom_line() +
  ggtitle(paste(c2, "cycle", c1, "EC50", sep = " ")) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x=element_blank()) +
  scale_color_manual(values = colors)

c2 <- "Light"

p5 <- data.frame(stage, SPF_light_RQ, GF_light_RQ) %>% melt() %>% 
  ggplot(aes(x=stage, y=value, group=variable, color=variable)) +
  geom_line() +
  ggtitle(paste(c2, "cycle", c1, "EC50", sep = " ")) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x=element_blank()) +
  scale_color_manual(values = colors)

c2 <- "Dark"

p6 <- data.frame(stage, SPF_dark_RQ, GF_dark_RQ) %>% melt() %>% 
  ggplot(aes(x=stage, y=value, group=variable, color=variable)) +
  geom_line() +
  ggtitle(paste(c2, "cycle", c1, "EC50", sep = " ")) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x=element_blank()) +
  scale_color_manual(values = colors)

ggarrange(p1,p2,p3,p4,p5,p6, common.legend = T) %>% 
  annotate_figure(left = text_grob("value"))

VO2

# VO2 ####

cuts <- seq(0,nrow(total_run),240)
cuts <- sapply(cuts, function(x) x+1)

channel = "VO2"

p <- graph_groups(total_run, channel, title = paste0(channel, " before and after adding HF diet"))
p + geom_rect(aes(xmin=cuts[6], xmax=cuts[7], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_rect(aes(xmin=cuts[8], xmax=cuts[9], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_rect(aes(xmin=cuts[10], xmax=cuts[11], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_rect(aes(xmin=cuts[12], xmax=cuts[13], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_rect(aes(xmin=cuts[14], xmax=cuts[15], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_line(size=1) +
  geom_vline(xintercept = 1500)

# p1 <- graph_groups(total2_pt1, channel, title = paste0(channel, " before HF diet"))
# p2 <- graph_groups(total2_pt2, channel, title = paste0(channel, " immediately after HF diet"))
# p3 <- graph_groups(total2_pt3, channel, title = paste0(channel, " 1 month after HF diet"))
p1 <- graph_groups(total2_pt1, channel)
p2 <- graph_groups(total2_pt2, channel)
p3 <- graph_groups(total2_pt3, channel)
leg <- get_legend(p1)
p1 <- p1 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
p2 <- p2 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
p3 <- p3 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
f <- ggarrange(p1, p2, p3, leg)
annotate_figure(f, top = text_grob(paste0(channel, ", before, after, and a month after adding HF diet"), face = "bold"), left = text_grob(channel), bottom = text_grob("observation"))

Before HF diet

# pt1 VO2 ####

pt = "before HF diet"
cycle = "Total"
total_pt1.VO2.rcf <- rcf_pipe(total_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.8746158 4.9918380
## 
## $cut
## [1] 0.08234444

## [1] "SPF EC50: 1.75001617929315"
## [1] "GF EC50: 1.60877208757266"
t <- total_pt1.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.1_VO2_M_6 = 0.42252, p-value = 0.2344
## alternative hypothesis: highest value 2.08865368747895 is an outlier
dixon.test(t[[1]][-which.max(t[[1]])])
## 
##  Dixon test for outliers
## 
## data:  t[[1]][-which.max(t[[1]])]
## Q.cohort1.1_VO2_M_9 = 0.43748, p-value = 0.2774
## alternative hypothesis: lowest value 1.46653571141594 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.1_VO2_M_10 = 0.049739, p-value = 0.2924
## alternative hypothesis: lowest value 1.42989814868382 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.97779, p-value = 0.952
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.91597, p-value = 0.3599
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 2.3404, num df = 8, denom df = 8, p-value = 0.2504
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.5279073 10.3753815
## sample estimates:
## ratio of variances 
##            2.34035
# t.test(t[[1]], t[[2]])
# If var.test is not significant, set var.equal to true!
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 2.0526, df = 16, p-value = 0.05684
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.00474476  0.29423557
## sample estimates:
## mean of x mean of y 
##  1.749973  1.605228
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7152 -1.5104 -0.0060  0.7461  5.6158 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  27.9665     5.7561   4.859 0.000174 ***
## df[[2]]       0.8346     3.4159   0.244 0.810083    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.297 on 16 degrees of freedom
## Multiple R-squared:  0.003717,   Adjusted R-squared:  -0.05855 
## F-statistic: 0.0597 on 1 and 16 DF,  p-value: 0.8101
## 
## y = 28 + 0.83x, r^2 = 0.00372

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.42989815   1.75896428   0.32906613 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  14.44704950   1.63768933   1.60522772   0.03858337   0.08897341   0.01339809 
##      std.dev     coef.var 
##   0.11575011   0.07210822 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.46653571   2.08865369   0.62211798 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  15.74975817   1.75334199   1.74997313   0.05902563   0.13611334   0.03135622 
##      std.dev     coef.var 
##   0.17707688   0.10118834 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 276.00000000  30.00000000  30.66666667   0.68211273   1.57295478   4.18750000 
##      std.dev     coef.var 
##   2.04633819   0.06672842 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  25.60000000  30.40000000   4.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 252.60000000  27.80000000  28.06666667   0.53800041   1.24063118   2.60500000 
##      std.dev     coef.var 
##   1.61400124   0.05750598 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.8135 0.3805
##       16               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0943 0.09428   4.213 0.0568 .
## Residuals   16 0.3580 0.02238                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  30.42  30.420   8.957 0.00861 **
## Residuals   16  54.34   3.396                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## df[[1]]      1 0.0004462 0.0004462   20.34 0.000356 ***
## Residuals   16 0.0003509 0.0000219                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.09428 0.09428   5.067 0.0398 *
## df[[3]]      1 0.07893 0.07893   4.242 0.0572 .
## Residuals   15 0.27910 0.01861                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value   Pr(>F)   
## (Intercept) 0.010936  1  0.5878 0.455182   
## df[[1]]     0.171531  1  9.2187 0.008334 **
## df[[3]]     0.078931  1  4.2421 0.057229 . 
## Residuals   0.279103 15                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.00833428507313233"
cycle = "Light"
light_pt1.VO2.rcf <- rcf_pipe(light_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.8746158 2.9373260
## 
## $cut
## [1] 0.0412542

## [1] "SPF EC50: 1.58906026308396"
## [1] "GF EC50: 1.44574161742324"
t <- light_pt1.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.1_VO2_M_6 = 0.19725, p-value = 0.9629
## alternative hypothesis: highest value 1.79542307671054 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort1.1_VO2_M_4 = 0.076364, p-value = 0.4349
## alternative hypothesis: lowest value 1.26352543847839 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.92042, p-value = 0.3958
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.82751, p-value = 0.04187
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 2.499, num df = 8, denom df = 8, p-value = 0.2167
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.5636947 11.0787387
## sample estimates:
## ratio of variances 
##           2.499005
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 2.2009, df = 16, p-value = 0.04277
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.005185257 0.276681933
## sample estimates:
## mean of x mean of y 
##  1.568389  1.427456
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7341 -1.5320 -0.0815  0.8395  5.6255 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   28.913      5.584   5.177 9.17e-05 ***
## df[[2]]        0.303      3.711   0.082    0.936    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.301 on 16 degrees of freedom
## Multiple R-squared:  0.0004167,  Adjusted R-squared:  -0.06206 
## F-statistic: 0.00667 on 1 and 16 DF,  p-value: 0.9359
## 
## y = 29 + 0.3x, r^2 = 0.000417

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.26352544   1.52389696   0.26037152 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  12.84710286   1.48944874   1.42745587   0.03423306   0.07894159   0.01054712 
##      std.dev     coef.var 
##   0.10269919   0.07194562 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.34191921   1.79542308   0.45350387 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  14.11550522   1.63088754   1.56838947   0.05411646   0.12479277   0.02635732 
##      std.dev     coef.var 
##   0.16234937   0.10351343 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 276.00000000  30.00000000  30.66666667   0.68211273   1.57295478   4.18750000 
##      std.dev     coef.var 
##   2.04633819   0.06672842 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  25.60000000  30.40000000   4.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 252.60000000  27.80000000  28.06666667   0.53800041   1.24063118   2.60500000 
##      std.dev     coef.var 
##   1.61400124   0.05750598 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   1.482 0.2411
##       16               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.08938 0.08938   4.844 0.0428 *
## Residuals   16 0.29524 0.01845                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  30.42  30.420   8.957 0.00861 **
## Residuals   16  54.34   3.396                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## df[[1]]      1 0.0003909 0.0003909   19.97 0.000387 ***
## Residuals   16 0.0003131 0.0000196                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.08938 0.08938   5.636 0.0314 *
## df[[3]]      1 0.05736 0.05736   3.617 0.0766 .
## Residuals   15 0.23788 0.01586                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value   Pr(>F)   
## (Intercept) 0.010671  1  0.6729 0.424895   
## df[[1]]     0.146579  1  9.2430 0.008267 **
## df[[3]]     0.057359  1  3.6169 0.076576 . 
## Residuals   0.237876 15                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.00826660982644304"
cycle = "Dark"
dark_pt1.VO2.rcf <- rcf_pipe(dark_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 1.000101 4.991838
## 
## $cut
## [1] 0.07983474

## [1] "SPF EC50: 1.9341144311299"
## [1] "GF EC50: 1.83403636087834"
t <- dark_pt1.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.1_VO2_M_6 = 0.47365, p-value = 0.1482
## alternative hypothesis: highest value 2.37415711632301 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.1_VO2_M_10 = 0.097525, p-value = 0.5498
## alternative hypothesis: lowest value 1.59851055132634 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.94522, p-value = 0.6378
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.93588, p-value = 0.5393
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 2.1262, num df = 8, denom df = 8, p-value = 0.3065
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4796022 9.4260035
## sample estimates:
## ratio of variances 
##           2.126201
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 1.6297, df = 16, p-value = 0.1227
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.04176041  0.31941330
## sample estimates:
## mean of x mean of y 
##  1.953008  1.814181
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6744 -1.3941  0.0982  0.7798  5.6352 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   26.732      5.541   4.824 0.000187 ***
## df[[2]]        1.399      2.928   0.478 0.639377    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.285 on 16 degrees of freedom
## Multiple R-squared:  0.01406,    Adjusted R-squared:  -0.04756 
## F-statistic: 0.2281 on 1 and 16 DF,  p-value: 0.6394
## 
## y = 27 + 1.4x, r^2 = 0.0141

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.59851055   2.00713925   0.40862869 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  16.32762964   1.82944446   1.81418107   0.04817938   0.11110185   0.02089127 
##      std.dev     coef.var 
##   0.14453814   0.07967129 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.60328120   2.37415712   0.77087591 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  17.57706766   1.93427500   1.95300752   0.07025276   0.16200316   0.04441905 
##      std.dev     coef.var 
##   0.21075829   0.10791473 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 276.00000000  30.00000000  30.66666667   0.68211273   1.57295478   4.18750000 
##      std.dev     coef.var 
##   2.04633819   0.06672842 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  25.60000000  30.40000000   4.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 252.60000000  27.80000000  28.06666667   0.53800041   1.24063118   2.60500000 
##      std.dev     coef.var 
##   1.61400124   0.05750598 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.2463 0.6264
##       16               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0867 0.08673   2.656  0.123
## Residuals   16 0.5225 0.03266               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  30.42  30.420   8.957 0.00861 **
## Residuals   16  54.34   3.396                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.0004835 0.0004835   15.65 0.00113 **
## Residuals   16 0.0004944 0.0000309                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0867 0.08673   3.176 0.0950 .
## df[[3]]      1 0.1128 0.11284   4.132 0.0602 .
## Residuals   15 0.4096 0.02731                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.00997  1  0.3651 0.55474  
## df[[1]]     0.19101  1  6.9941 0.01839 *
## df[[3]]     0.11284  1  4.1320 0.06017 .
## Residuals   0.40964 15                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0183871212047179"

Immediately after HF diet

# pt2 VO2 ####

# delete these once I add cohort1_2!!
covar_SPF_pt1 <- c(cohort2_1_covar[[1]][c(6,7,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
covar_GF_pt1 <- c(cohort2_1_covar[[1]][c(1,2,3,4,5)], cohort3_1_covar[[1]][c(1,2,3)])

pt = "immediately after HF diet"
cycle = "Total"
total_pt2.VO2.rcf <- rcf_pipe(total_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.8753851 3.3065780
## 
## $cut
## [1] 0.04862386

## [1] "SPF EC50: 1.81300642705308"
## [1] "GF EC50: 1.6968627087087"
t <- total_pt2.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.2_VO2_M_6 = 0.39782, p-value = 0.2693
## alternative hypothesis: highest value 2.33215178662138 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort2.2_VO2_M_1 = 0.39801, p-value = 0.3661
## alternative hypothesis: lowest value 1.46610566358194 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.94242, p-value = 0.6606
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.94787, p-value = 0.6897
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 5.7567, num df = 6, denom df = 7, p-value = 0.03675
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   1.124668 32.787226
## sample estimates:
## ratio of variances 
##           5.756719
t.test(t[[1]], t[[2]], var.equal = FALSE) # F-test p-val: 0.03675
## 
##  Welch Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 1.2272, df = 7.808, p-value = 0.2555
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1244126  0.4048969
## sample estimates:
## mean of x mean of y 
##  1.829071  1.688829
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6643 -1.2707 -0.2855  0.5268  5.1774 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   27.725      4.901   5.657 7.84e-05 ***
## df[[2]]        1.141      2.774   0.411    0.688    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.23 on 13 degrees of freedom
## Multiple R-squared:  0.01285,    Adjusted R-squared:  -0.06309 
## F-statistic: 0.1692 on 1 and 13 DF,  p-value: 0.6875
## 
## y = 28 + 1.1x, r^2 = 0.0128

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000   1.46610566   1.83829826   0.37219260 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  13.51063076   1.70811396   1.68882884   0.04150956   0.09815452   0.01378435 
##      std.dev     coef.var 
##   0.11740677   0.06951964 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   7.00000000   0.00000000   0.00000000   1.40283486   2.33215179   0.92931693 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  12.80349695   1.79453176   1.82907099   0.10647108   0.26052535   0.07935264 
##      std.dev     coef.var 
##   0.28169600   0.15401043 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 246.50000000  30.00000000  30.81250000   0.75555975   1.78661491   4.56696429 
##      std.dev     coef.var 
##   2.13704569   0.06935645 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##    7.0000000    0.0000000    0.0000000   26.3000000   30.4000000    4.1000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  199.4000000   28.1000000   28.4857143    0.5629127    1.3773978    2.2180952 
##      std.dev     coef.var 
##    1.4893271    0.0522833 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.8247 0.1998
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0734 0.07343   1.667  0.219
## Residuals   13 0.5726 0.04405               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1  20.21  20.212   5.803 0.0315 *
## Residuals   13  45.28   3.483                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0003232 0.0003232   7.064 0.0197 *
## Residuals   13 0.0005949 0.0000458                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0734 0.07343   1.805  0.204
## df[[3]]      1 0.0845 0.08446   2.076  0.175
## Residuals   12 0.4881 0.04068               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.00608  1  0.1494 0.70590  
## df[[1]]     0.14959  1  3.6773 0.07927 .
## df[[3]]     0.08446  1  2.0763 0.17518  
## Residuals   0.48814 12                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)


cycle = "Light"
light_pt2.VO2.rcf <- rcf_pipe(light_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.1013106 3.2033890
## 
## $cut
## [1] 0.06204157

## [1] "SPF EC50: 1.71309371902381"
## [1] "GF EC50: 1.59241129302527"
t <- light_pt2.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.2_VO2_M_6 = 0.40248, p-value = 0.2596
## alternative hypothesis: highest value 2.05964348633281 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.2_VO2_M_10 = 0.25228, p-value = 0.8294
## alternative hypothesis: lowest value 1.39075334998537 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.95817, p-value = 0.8029
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.91587, p-value = 0.3973
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 2.8859, num df = 6, denom df = 7, p-value = 0.1917
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.5638061 16.4365354
## sample estimates:
## ratio of variances 
##           2.885896
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 1.5625, df = 13, p-value = 0.1422
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.04971928  0.30961818
## sample estimates:
## mean of x mean of y 
##  1.713522  1.583573
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5318 -1.2153 -0.2414  0.5103  5.2058 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   28.068      5.855   4.794 0.000351 ***
## df[[2]]        1.009      3.543   0.285 0.780371    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.238 on 13 degrees of freedom
## Multiple R-squared:  0.006196,   Adjusted R-squared:  -0.07025 
## F-statistic: 0.08104 on 1 and 13 DF,  p-value: 0.7804
## 
## y = 28 + 1x, r^2 = 0.0062

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000   1.39075335   1.71114273   0.32038938 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  12.66858433   1.60455785   1.58357304   0.04154122   0.09822937   0.01380538 
##      std.dev     coef.var 
##   0.11749630   0.07419696 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   7.00000000   0.00000000   0.00000000   1.40092810   2.05964349   0.65871539 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  11.99465742   1.72520434   1.71352249   0.07544240   0.18460091   0.03984089 
##      std.dev     coef.var 
##   0.19960184   0.11648627 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 246.50000000  30.00000000  30.81250000   0.75555975   1.78661491   4.56696429 
##      std.dev     coef.var 
##   2.13704569   0.06935645 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##    7.0000000    0.0000000    0.0000000   26.3000000   30.4000000    4.1000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  199.4000000   28.1000000   28.4857143    0.5629127    1.3773978    2.2180952 
##      std.dev     coef.var 
##    1.4893271    0.0522833 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.4183  0.529
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0630 0.06304   2.442  0.142
## Residuals   13 0.3357 0.02582               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1  20.21  20.212   5.803 0.0315 *
## Residuals   13  45.28   3.483                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.0002788 2.788e-04   10.46 0.00653 **
## Residuals   13 0.0003466 2.666e-05                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df  Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.06304 0.06304   2.665  0.129
## df[[3]]      1 0.05177 0.05177   2.188  0.165
## Residuals   12 0.28391 0.02366               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.013909  1  0.5879 0.45806  
## df[[1]]     0.112346  1  4.7485 0.04997 *
## df[[3]]     0.051772  1  2.1882 0.16483  
## Residuals   0.283911 12                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0499736001563136"
cycle = "Dark"
dark_pt2.VO2.rcf <- rcf_pipe(dark_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 1.032885 3.306578
## 
## $cut
## [1] 0.04547386

## [1] "SPF EC50: 1.92194220183729"
## [1] "GF EC50: 1.84062090652473"
t <- dark_pt2.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.2_VO2_M_6 = 0.31187, p-value = 0.4956
## alternative hypothesis: highest value 2.60800348939703 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.2_VO2_M_9 = 0.37516, p-value = 0.4268
## alternative hypothesis: lowest value 1.46852607705669 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.97884, p-value = 0.9537
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.84905, p-value = 0.09319
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 4.5986, num df = 6, denom df = 7, p-value = 0.06562
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.898416 26.191357
## sample estimates:
## ratio of variances 
##           4.598629
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 0.93562, df = 13, p-value = 0.3665
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.1880892  0.4754632
## sample estimates:
## mean of x mean of y 
##  1.955075  1.811388
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6876 -1.2898 -0.2909  0.5627  5.2043 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.3373     3.8384   7.383 5.33e-06 ***
## df[[2]]       0.7396     2.0202   0.366     0.72    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.233 on 13 degrees of freedom
## Multiple R-squared:  0.01021,    Adjusted R-squared:  -0.06593 
## F-statistic: 0.134 on 1 and 13 DF,  p-value: 0.7202
## 
## y = 28 + 0.74x, r^2 = 0.0102

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000   1.46852608   1.97171082   0.50318475 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  14.49110787   1.89549878   1.81138848   0.06431407   0.15207860   0.03309039 
##      std.dev     coef.var 
##   0.18190765   0.10042443 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##    7.0000000    0.0000000    0.0000000    1.3997469    2.6080035    1.2082566 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##   13.6855283    1.8564543    1.9550755    0.1474403    0.3607734    0.1521705 
##      std.dev     coef.var 
##    0.3900903    0.1995270 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 246.50000000  30.00000000  30.81250000   0.75555975   1.78661491   4.56696429 
##      std.dev     coef.var 
##   2.13704569   0.06935645 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##    7.0000000    0.0000000    0.0000000   26.3000000   30.4000000    4.1000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  199.4000000   28.1000000   28.4857143    0.5629127    1.3773978    2.2180952 
##      std.dev     coef.var 
##    1.4893271    0.0522833 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.9939 0.1814
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0771 0.07708   0.875  0.367
## Residuals   13 1.1447 0.08805               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1  20.21  20.212   5.803 0.0315 *
## Residuals   13  45.28   3.483                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0003571 0.0003571   3.648 0.0785 .
## Residuals   13 0.0012727 0.0000979                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0771 0.07708   0.887  0.365
## df[[3]]      1 0.1023 0.10226   1.177  0.299
## Residuals   12 1.0424 0.08687               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value Pr(>F)
## (Intercept) 0.00571  1  0.0657 0.8020
## df[[1]]     0.16687  1  1.9210 0.1910
## df[[3]]     0.10226  1  1.1773 0.2992
## Residuals   1.04239 12
check_ANCOVA(res)

1 month after HF diet

# pt3 VO2 ####

# cohort2.3_VO2_M_1 is an outlier.
total_pt3 %<>% select(-cohort2.3_VO2_M_1)
light_pt3 %<>% select(-cohort2.3_VO2_M_1)
dark_pt3 %<>% select(-cohort2.3_VO2_M_1)
covar_SPF_pt3 <- c(cohort1_3_covar[[1]][c(2,3)], cohort2_3_covar[[1]][c(4,5,6,7)], cohort3_3_covar[[1]][c(4,5,6)])
covar_GF_pt3 <- c(cohort1_3_covar[[1]][1], cohort2_3_covar[[1]][c(2,3)], cohort3_3_covar[[1]][c(1,2,3)])

pt = "1 month after HF diet"
cycle = "Total"
total_pt3.VO2.rcf <- rcf_pipe(total_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.1282987 3.6305530
## 
## $cut
## [1] 0.07004509

## [1] "SPF EC50: 1.86394965004008"
## [1] "GF EC50: 1.68055953099069"
t <- total_pt3.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.3_VO2_M_9 = 0.13836, p-value = 0.7538
## alternative hypothesis: lowest value 1.64628551511344 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.3_VO2_M_3 = 0.31403, p-value = 0.6131
## alternative hypothesis: highest value 1.7862408128281 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.95564, p-value = 0.7518
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.99256, p-value = 0.9945
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 6.0101, num df = 8, denom df = 5, p-value = 0.06389
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.889446 28.952498
## sample estimates:
## ratio of variances 
##           6.010139
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 2.758, df = 13, p-value = 0.01629
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.0401974 0.3308189
## sample estimates:
## mean of x mean of y 
##  1.875537  1.690029
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7604 -2.7232 -0.9105  2.5311  5.9239 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   30.707      9.858   3.115  0.00821 **
## df[[2]]        1.014      5.454   0.186  0.85538   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.16 on 13 degrees of freedom
## Multiple R-squared:  0.002652,   Adjusted R-squared:  -0.07407 
## F-statistic: 0.03457 on 1 and 13 DF,  p-value: 0.8554
## 
## y = 31 + 1x, r^2 = 0.00265

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  6.000000000  0.000000000  0.000000000  1.605187543  1.786240813  0.181053270 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 10.140175624  1.683826709  1.690029271  0.025783815  0.066279407  0.003988831 
##      std.dev     coef.var 
##  0.063157191  0.037370472 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.64628552   2.09318823   0.44690271 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  16.87983684   1.87746493   1.87553743   0.05161118   0.11901560   0.02397343 
##      std.dev     coef.var 
##   0.15483355   0.08255423 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   6.00000000   0.00000000   0.00000000  28.70000000  33.20000000   4.50000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 179.60000000  29.35000000  29.93333333   0.67806915   1.74303225   2.75866667 
##      std.dev     coef.var 
##   1.66092344   0.05548742 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  30.90000000  38.30000000   7.40000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 308.40000000  34.90000000  34.26666667   0.82276634   1.89730257   6.09250000 
##      std.dev     coef.var 
##   2.46829901   0.07203207 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.6773 0.07739 .
##       13                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.1239 0.12389   7.607 0.0163 *
## Residuals   13 0.2117 0.01629                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  67.60   67.60   14.05 0.00243 **
## Residuals   13  62.53    4.81                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.0000070 7.030e-06   0.151  0.704
## Residuals   13 0.0006059 4.661e-05               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.1239 0.12389   13.84 0.00293 **
## df[[3]]      1 0.1043 0.10428   11.65 0.00515 **
## Residuals   12 0.1075 0.00895                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.58516  1  65.349 3.379e-06 ***
## df[[1]]     0.22728  1  25.382 0.0002903 ***
## df[[3]]     0.10428  1  11.646 0.0051489 ** 
## Residuals   0.10745 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV and Covar,"
cycle = "Light"
light_pt3.VO2.rcf <- rcf_pipe(light_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.1282987 3.2196870
## 
## $cut
## [1] 0.06182777

## [1] "SPF EC50: 1.69008429724231"
## [1] "GF EC50: 1.45382056067872"
t <- light_pt3.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.3_VO2_M_9 = 0.20485, p-value = 0.9297
## alternative hypothesis: lowest value 1.43304699918151 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort2.3_VO2_M_3 = 0.11382, p-value = 0.5661
## alternative hypothesis: highest value 1.53061348786286 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.97217, p-value = 0.9126
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.93438, p-value = 0.6144
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 6.0985, num df = 8, denom df = 5, p-value = 0.06199
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.9025158 29.3779346
## sample estimates:
## ratio of variances 
##           6.098454
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 3.5555, df = 13, p-value = 0.003519
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.08659019 0.35475341
## sample estimates:
## mean of x mean of y 
##  1.677844  1.457173
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.666 -2.449 -1.202  2.489  6.212 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   28.012      8.379   3.343  0.00529 **
## df[[2]]        2.844      5.247   0.542  0.59695   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.129 on 13 degrees of freedom
## Multiple R-squared:  0.0221, Adjusted R-squared:  -0.05312 
## F-statistic: 0.2938 on 1 and 13 DF,  p-value: 0.597
## 
## y = 28 + 2.8x, r^2 = 0.0221

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  6.000000000  0.000000000  0.000000000  1.387532895  1.530613488  0.143080593 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  8.743035582  1.450497011  1.457172597  0.023634555  0.060754559  0.003351553 
##      std.dev     coef.var 
##  0.057892601  0.039729405 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.43304700   1.89381736   0.46077036 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  15.10059959   1.71904878   1.67784440   0.04765535   0.10989344   0.02043929 
##      std.dev     coef.var 
##   0.14296606   0.08520818 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   6.00000000   0.00000000   0.00000000  28.70000000  33.20000000   4.50000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 179.60000000  29.35000000  29.93333333   0.67806915   1.74303225   2.75866667 
##      std.dev     coef.var 
##   1.66092344   0.05548742 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  30.90000000  38.30000000   7.40000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 308.40000000  34.90000000  34.26666667   0.82276634   1.89730257   6.09250000 
##      std.dev     coef.var 
##   2.46829901   0.07203207 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.9704 0.1838
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.1753 0.17531   12.64 0.00352 **
## Residuals   13 0.1803 0.01387                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  67.60   67.60   14.05 0.00243 **
## Residuals   13  62.53    4.81                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.0000011 1.100e-06   0.029  0.868
## Residuals   13 0.0004959 3.814e-05               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df  Sum Sq Mean Sq F value   Pr(>F)    
## df[[1]]      1 0.17531 0.17531   24.53 0.000335 ***
## df[[3]]      1 0.09452 0.09452   13.23 0.003407 ** 
## Residuals   12 0.08575 0.00715                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.47390  1  66.318 3.132e-06 ***
## df[[1]]     0.26197  1  36.660 5.717e-05 ***
## df[[3]]     0.09452  1  13.227  0.003407 ** 
## Residuals   0.08575 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV and Covar,"
cycle = "Dark"
dark_pt3.VO2.rcf <- rcf_pipe(dark_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 1.133380 3.630553
## 
## $cut
## [1] 0.04994346

## [1] "SPF EC50: 2.05497861966052"
## [1] "GF EC50: 1.93321699617181"
t <- dark_pt3.VO2.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.3_VO2_M_12 = 0.00219, p-value = 0.01876
## alternative hypothesis: highest value 2.35630507609501 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort1.3_VO2_M_10 = 0.2955, p-value = 0.6757
## alternative hypothesis: lowest value 1.76718655653742 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.90042, p-value = 0.2544
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.91308, p-value = 0.457
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 3.0977, num df = 8, denom df = 5, p-value = 0.2285
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.4584277 14.9223546
## sample estimates:
## ratio of variances 
##           3.097675
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 2.1269, df = 13, p-value = 0.05315
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.002649689  0.339813505
## sample estimates:
## mean of x mean of y 
##  2.088823  1.920241
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7689 -2.7219 -0.8649  2.5665  5.8477 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  31.2124    10.1851   3.064  0.00904 **
## df[[2]]       0.6535     5.0225   0.130  0.89847   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.162 on 13 degrees of freedom
## Multiple R-squared:  0.001301,   Adjusted R-squared:  -0.07552 
## F-statistic: 0.01693 on 1 and 13 DF,  p-value: 0.8985
## 
## y = 31 + 0.65x, r^2 = 0.0013

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  6.000000000  0.000000000  0.000000000  1.767186557  2.016871555  0.249684998 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 11.521444188  1.941350371  1.920240698  0.040563066  0.104270682  0.009872174 
##      std.dev     coef.var 
##  0.099358815  0.051742896 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000   1.87848131   2.35630508   0.47782377 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  18.79940346   2.06479867   2.08882261   0.05829121   0.13441978   0.03058079 
##      std.dev     coef.var 
##   0.17487364   0.08371876 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   6.00000000   0.00000000   0.00000000  28.70000000  33.20000000   4.50000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 179.60000000  29.35000000  29.93333333   0.67806915   1.74303225   2.75866667 
##      std.dev     coef.var 
##   1.66092344   0.05548742 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  30.90000000  38.30000000   7.40000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 308.40000000  34.90000000  34.26666667   0.82276634   1.89730257   6.09250000 
##      std.dev     coef.var 
##   2.46829901   0.07203207 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.0775 0.3182
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.1023 0.10231   4.524 0.0531 .
## Residuals   13 0.2940 0.02262                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  67.60   67.60   14.05 0.00243 **
## Residuals   13  62.53    4.81                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.0000271 2.706e-05   0.471  0.504
## Residuals   13 0.0007465 5.742e-05               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df  Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.10231 0.10231   6.015 0.0305 *
## df[[3]]      1 0.08989 0.08989   5.285 0.0403 *
## Residuals   12 0.20412 0.01701                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.64393  1 37.8564 4.926e-05 ***
## df[[1]]     0.19169  1 11.2691  0.005706 ** 
## df[[3]]     0.08989  1  5.2846  0.040283 *  
## Residuals   0.20412 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV and Covar,"

RQ

# RQ ####

channel = "RQ"

total_run %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>% 
  mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.5, NA, .))
total2_pt2 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>% 
  mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.1, NA, .))
total2_pt3 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .))
total_pt2 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>% 
  mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.1, NA, .))
total_pt3 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .))
light_pt2 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>% 
  mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.1, NA, .))
dark_pt2 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .)) %>% 
  mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. > 1.1, NA, .))
light_pt3 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .))
dark_pt3 %<>% mutate_at(dplyr::vars(contains("RQ")), ~ ifelse(. < 0.5, NA, .))

p <- graph_groups(total_run, channel, title = paste0(channel, " before and after adding HF diet"))
p + geom_rect(aes(xmin=cuts[6], xmax=cuts[7], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_rect(aes(xmin=cuts[8], xmax=cuts[9], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_rect(aes(xmin=cuts[10], xmax=cuts[11], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_rect(aes(xmin=cuts[12], xmax=cuts[13], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_rect(aes(xmin=cuts[14], xmax=cuts[15], ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
  geom_line(size=1) +
  geom_vline(xintercept = 1500)

# p1 <- graph_groups(total2_pt1, channel, title = paste0(channel, " before HF diet"))
# p2 <- graph_groups(total2_pt2, channel, title = paste0(channel, " immediately after HF diet"))
# p3 <- graph_groups(total2_pt3, channel, title = paste0(channel, " 1 month after HF diet"))
p1 <- graph_groups(total2_pt1, channel)
p2 <- graph_groups(total2_pt2, channel)
p3 <- graph_groups(total2_pt3, channel)
leg <- get_legend(p1)
p1 <- p1 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
p2 <- p2 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
p3 <- p3 + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank())
f <- ggarrange(p1, p2, p3, leg)
annotate_figure(f, top = text_grob(paste0(channel, ", before, after, and a month after adding HF diet"), face = "bold"), left = text_grob(channel), bottom = text_grob("observation"))

Before HF diet

# pt1 RQ ####

covar_SPF_pt1 <- c(cohort1_1_covar[[1]][c(2,3)], cohort2_1_covar[[1]][c(6,7,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
covar_GF_pt1 <- c(cohort1_1_covar[[1]][1], cohort2_1_covar[[1]][c(1,2,3,4,5)], cohort3_1_covar[[1]][c(1,2,3)])

pt = "before HF diet"
cycle = "Total"
total_pt1.RQ.rcf <- rcf_pipe(total_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.2911567 1.3067190
## 
## $cut
## [1] 0.02031125

## [1] "SPF EC50: 0.86888794275416"
## [1] "GF EC50: 0.826388799709833"
t <- total_pt1.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.1_RQ_M_9 = 0.48246, p-value = 0.1359
## alternative hypothesis: lowest value 0.823158188054464 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.1_RQ_M_9 = 0.56455, p-value = 0.0536
## alternative hypothesis: lowest value 0.760142288370411 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.95552, p-value = 0.7505
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.86605, p-value = 0.1115
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 0.73521, num df = 8, denom df = 8, p-value = 0.6738
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1658385 3.2593552
## sample estimates:
## ratio of variances 
##          0.7352051
# t.test(t[[1]], t[[2]])
# If var.test is not significant, set var.equal to true!
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 3.2744, df = 16, p-value = 0.00477
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.01572764 0.07348722
## sample estimates:
## mean of x mean of y 
## 0.8719954 0.8273880
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4906 -1.1167 -0.3626  1.0107  5.2180 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    51.85      11.83   4.382 0.000465 ***
## df[[2]]       -26.46      13.91  -1.902 0.075384 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.079 on 16 degrees of freedom
## Multiple R-squared:  0.1843, Adjusted R-squared:  0.1334 
## F-statistic: 3.616 on 1 and 16 DF,  p-value: 0.07538
## 
## y = 52 + -26x, r^2 = 0.184

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 9.0000000000 0.0000000000 0.0000000000 0.7601422884 0.8762337005 0.1160914122 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 7.4464918248 0.8317564691 0.8273879805 0.0103419439 0.0238485654 0.0009626022 
##      std.dev     coef.var 
## 0.0310258317 0.0374985284 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  9.000000000  0.000000000  0.000000000  0.823158188  0.914947835  0.091789646 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  7.847958723  0.864325085  0.871995414  0.008867607  0.020448738  0.000707710 
##      std.dev     coef.var 
##  0.026602820  0.030507982 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 276.00000000  30.00000000  30.66666667   0.68211273   1.57295478   4.18750000 
##      std.dev     coef.var 
##   2.04633819   0.06672842 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  25.60000000  30.40000000   4.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 252.60000000  27.80000000  28.06666667   0.53800041   1.24063118   2.60500000 
##      std.dev     coef.var 
##   1.61400124   0.05750598 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0132   0.91
##       16               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df   Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.008954 0.008954   10.72 0.00477 **
## Residuals   16 0.013362 0.000835                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  30.42  30.420   8.957 0.00861 **
## Residuals   16  54.34   3.396                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## df[[1]]      1 7.502e-05 7.502e-05   16.98 0.000801 ***
## Residuals   16 7.069e-05 4.420e-06                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq  Mean Sq F value Pr(>F)   
## df[[1]]      1 0.008954 0.008954  10.117 0.0062 **
## df[[3]]      1 0.000087 0.000087   0.098 0.7587   
## Residuals   15 0.013276 0.000885                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.043067  1 48.6603 4.458e-06 ***
## df[[1]]     0.004927  1  5.5667   0.03228 *  
## df[[3]]     0.000087  1  0.0979   0.75872    
## Residuals   0.013276 15                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0322840364003511"
cycle = "Light"
light_pt1.RQ.rcf <- rcf_pipe(light_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.4186665 0.9756889
## 
## $cut
## [1] 0.01114045

## [1] "SPF EC50: 0.808304086759372"
## [1] "GF EC50: 0.756862038448115"
t <- light_pt1.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.1_RQ_M_11 = 0.12043, p-value = 0.6669
## alternative hypothesis: highest value 0.855640198857886 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.1_RQ_M_9 = 0.45221, p-value = 0.181
## alternative hypothesis: lowest value 0.70776270316176 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.89747, p-value = 0.2377
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.91972, p-value = 0.3899
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 2.4598, num df = 8, denom df = 8, p-value = 0.2245
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.5548409 10.9047281
## sample estimates:
## ratio of variances 
##           2.459754
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 3.2633, df = 16, p-value = 0.004882
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.01710230 0.08051796
## sample estimates:
## mean of x mean of y 
## 0.8054247 0.7566145
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1339 -0.8929 -0.3644  0.7451  4.9434 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   47.952      9.955   4.817  0.00019 ***
## df[[2]]      -23.796     12.731  -1.869  0.08001 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.085 on 16 degrees of freedom
## Multiple R-squared:  0.1792, Adjusted R-squared:  0.1279 
## F-statistic: 3.494 on 1 and 16 DF,  p-value: 0.08001
## 
## y = 48 + -24x, r^2 = 0.179

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 9.0000000000 0.0000000000 0.0000000000 0.7077627032 0.7865111501 0.0787484469 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 6.8095307367 0.7534788906 0.7566145263 0.0080413220 0.0185433218 0.0005819657 
##      std.dev     coef.var 
## 0.0241239660 0.0318840904 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  9.000000000  0.000000000  0.000000000  0.756513163  0.855640199  0.099127036 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  7.248821878  0.814399686  0.805424653  0.012611690  0.029082609  0.001431492 
##      std.dev     coef.var 
##  0.037835069  0.046975306 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 276.00000000  30.00000000  30.66666667   0.68211273   1.57295478   4.18750000 
##      std.dev     coef.var 
##   2.04633819   0.06672842 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  25.60000000  30.40000000   4.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 252.60000000  27.80000000  28.06666667   0.53800041   1.24063118   2.60500000 
##      std.dev     coef.var 
##   1.61400124   0.05750598 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.5522 0.07776 .
##       16                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df  Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.01072 0.010721   10.65 0.00488 **
## Residuals   16 0.01611 0.001007                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  30.42  30.420   8.957 0.00861 **
## Residuals   16  54.34   3.396                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## df[[1]]      1 7.257e-05 7.257e-05   17.99 0.000622 ***
## Residuals   16 6.454e-05 4.030e-06                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.010721 0.010721  10.036 0.00637 **
## df[[3]]      1 0.000083 0.000083   0.078 0.78371   
## Residuals   15 0.016024 0.001068                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.036251  1 33.9339 3.342e-05 ***
## df[[1]]     0.005996  1  5.6126   0.03168 *  
## df[[3]]     0.000083  1  0.0781   0.78371    
## Residuals   0.016024 15                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0316788344413115"
cycle = "Dark"
dark_pt1.RQ.rcf <- rcf_pipe(dark_pt1, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.2911567 1.3067190
## 
## $cut
## [1] 0.02031125

## [1] "SPF EC50: 0.905887016006714"
## [1] "GF EC50: 0.888108777077216"
t <- dark_pt1.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort3.1_RQ_M_12 = 0.44552, p-value = 0.1922
## alternative hypothesis: lowest value 0.853686998799925 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.1_RQ_M_9 = 0.19108, p-value = 0.9903
## alternative hypothesis: lowest value 0.820799851235852 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.96529, p-value = 0.8518
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.90587, p-value = 0.288
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 0.68623, num df = 8, denom df = 8, p-value = 0.6068
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1547903 3.0422174
## sample estimates:
## ratio of variances 
##          0.6862258
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 1.8123, df = 16, p-value = 0.08875
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.004346843  0.055560099
## sample estimates:
## mean of x mean of y 
## 0.9066351 0.8810285
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6795 -1.2565 -0.2520  0.8756  5.8246 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    41.77      15.33   2.725    0.015 *
## df[[2]]       -13.87      17.14  -0.809    0.430  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.256 on 16 degrees of freedom
## Multiple R-squared:  0.03933,    Adjusted R-squared:  -0.02071 
## F-statistic: 0.6551 on 1 and 16 DF,  p-value: 0.4302
## 
## y = 42 + -14x, r^2 = 0.0393

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  9.000000000  0.000000000  0.000000000  0.820799851  0.919892794  0.099092943 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  7.929256641  0.890881776  0.881028516  0.010881101  0.025091863  0.001065585 
##      std.dev     coef.var 
##  0.032643302  0.037051357 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 9.0000000000 0.0000000000 0.0000000000 0.8536869988 0.9504519265 0.0967649277 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 8.1597162900 0.9084909024 0.9066351433 0.0090137676 0.0207857854 0.0007312321 
##      std.dev     coef.var 
## 0.0270413029 0.0298260034 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 276.00000000  30.00000000  30.66666667   0.68211273   1.57295478   4.18750000 
##      std.dev     coef.var 
##   2.04633819   0.06672842 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  25.60000000  30.40000000   4.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 252.60000000  27.80000000  28.06666667   0.53800041   1.24063118   2.60500000 
##      std.dev     coef.var 
##   1.61400124   0.05750598 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   0.196 0.6639
##       16               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.002951 0.0029506   3.284 0.0888 .
## Residuals   16 0.014375 0.0008984                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  30.42  30.420   8.957 0.00861 **
## Residuals   16  54.34   3.396                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)   
## df[[1]]      1 5.744e-05 5.744e-05    13.1 0.0023 **
## Residuals   16 7.016e-05 4.390e-06                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.002951 0.0029506   3.093  0.099 .
## df[[3]]      1 0.000065 0.0000646   0.068  0.798  
## Residuals   15 0.014310 0.0009540                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.041245  1 43.2339 8.805e-06 ***
## df[[1]]     0.002334  1  2.4464    0.1386    
## df[[3]]     0.000065  1  0.0678    0.7982    
## Residuals   0.014310 15                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)

Immediately after HF diet

# pt2 RQ ####

# delete these once I add cohort1_2!!
covar_SPF_pt1 <- c(cohort2_1_covar[[1]][c(6,7,8,9)], cohort3_1_covar[[1]][c(4,5,6)])
covar_GF_pt1 <- c(cohort2_1_covar[[1]][c(1,2,3,4,5)], cohort3_1_covar[[1]][c(1,2,3)])

pt = "immediately after HF diet"
cycle = "Total"
total_pt2.RQ.rcf <- rcf_pipe(total_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5205333 1.0818660
## 
## $cut
## [1] 0.01122665

## [1] "SPF EC50: 0.889997042895635"
## [1] "GF EC50: 0.813386849541187"
t <- total_pt2.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.2_RQ_M_9 = 0.4003, p-value = 0.2641
## alternative hypothesis: lowest value 0.841189809990053 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.2_RQ_M_9 = 0.20498, p-value = 0.9789
## alternative hypothesis: lowest value 0.724630874890402 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.90654, p-value = 0.3725
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.91694, p-value = 0.4055
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 0.25911, num df = 6, denom df = 7, p-value = 0.1205
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.0506214 1.4757565
## sample estimates:
## ratio of variances 
##          0.2591105
t.test(t[[1]], t[[2]], var.equal = FALSE) # F-test p-val: 0.03675
## 
##  Welch Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 4.0728, df = 10.668, p-value = 0.001961
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.0358323 0.1208002
## sample estimates:
## mean of x mean of y 
## 0.8851837 0.8068675
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5710 -1.2369 -0.4368  1.0528  4.4375 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   44.242      8.279   5.344 0.000133 ***
## df[[2]]      -17.210      9.796  -1.757 0.102472    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.018 on 13 degrees of freedom
## Multiple R-squared:  0.1919, Adjusted R-squared:  0.1297 
## F-statistic: 3.086 on 1 and 13 DF,  p-value: 0.1025
## 
## y = 44 + -17x, r^2 = 0.192

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  8.000000000  0.000000000  0.000000000  0.724630875  0.860256238  0.135625364 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  6.454939799  0.819557590  0.806867475  0.016890152  0.039938864  0.002282218 
##      std.dev     coef.var 
##  0.047772565  0.059207449 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 7.0000000000 0.0000000000 0.0000000000 0.8411898100 0.9090888065 0.0678989965 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 6.1962861453 0.8931365901 0.8851837350 0.0091911972 0.0224900493 0.0005913467 
##      std.dev     coef.var 
## 0.0243176220 0.0274718355 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 246.50000000  30.00000000  30.81250000   0.75555975   1.78661491   4.56696429 
##      std.dev     coef.var 
##   2.13704569   0.06935645 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##    7.0000000    0.0000000    0.0000000   26.3000000   30.4000000    4.1000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  199.4000000   28.1000000   28.4857143    0.5629127    1.3773978    2.2180952 
##      std.dev     coef.var 
##    1.4893271    0.0522833 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.4961  0.243
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df  Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.02290 0.022898   15.25 0.00181 **
## Residuals   13 0.01952 0.001502                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1  20.21  20.212   5.803 0.0315 *
## Residuals   13  45.28   3.483                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## df[[1]]      1 8.775e-05 8.775e-05   18.31 0.000897 ***
## Residuals   13 6.230e-05 4.790e-06                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.022898 0.022898  14.114 0.00274 **
## df[[3]]      1 0.000055 0.000055   0.034 0.85740   
## Residuals   12 0.019469 0.001622                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.033509  1 20.6538 0.0006724 ***
## df[[1]]     0.014814  1  9.1309 0.0106284 *  
## df[[3]]     0.000055  1  0.0337 0.8573954    
## Residuals   0.019469 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0106283891377386"
cycle = "Light"
light_pt2.RQ.rcf <- rcf_pipe(light_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5205333 1.0163130
## 
## $cut
## [1] 0.009915594

## [1] "SPF EC50: 0.866984758589404"
## [1] "GF EC50: 0.787869144945168"
t <- light_pt2.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort3.2_RQ_M_12 = 0.31499, p-value = 0.4858
## alternative hypothesis: lowest value 0.810903257955022 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.2_RQ_M_9 = 0.39154, p-value = 0.3827
## alternative hypothesis: lowest value 0.70694525860472 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.95353, p-value = 0.7617
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.9858, p-value = 0.9858
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 0.5345, num df = 6, denom df = 7, p-value = 0.463
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1044234 3.0442349
## sample estimates:
## ratio of variances 
##           0.534501
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 3.5647, df = 13, p-value = 0.003457
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.03017021 0.12299306
## sample estimates:
## mean of x mean of y 
## 0.8638716 0.7872899
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5778 -1.2277 -0.5805  1.1171  4.2483 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   45.344      7.651   5.927 5.01e-05 ***
## df[[2]]      -18.976      9.275  -2.046   0.0616 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.952 on 13 degrees of freedom
## Multiple R-squared:  0.2435, Adjusted R-squared:  0.1854 
## F-statistic: 4.185 on 1 and 13 DF,  p-value: 0.06156
## 
## y = 45 + -19x, r^2 = 0.244

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  8.000000000  0.000000000  0.000000000  0.706945259  0.860038846  0.153093587 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  6.298319331  0.787148223  0.787289916  0.016562386  0.039163820  0.002194501 
##      std.dev     coef.var 
##  0.046845502  0.059502225 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  7.000000000  0.000000000  0.000000000  0.810903258  0.905564075  0.094660817 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  6.047100861  0.867761027  0.863871552  0.012944735  0.031674625  0.001172963 
##      std.dev     coef.var 
##  0.034248549  0.039645418 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 246.50000000  30.00000000  30.81250000   0.75555975   1.78661491   4.56696429 
##      std.dev     coef.var 
##   2.13704569   0.06935645 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##    7.0000000    0.0000000    0.0000000   26.3000000   30.4000000    4.1000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  199.4000000   28.1000000   28.4857143    0.5629127    1.3773978    2.2180952 
##      std.dev     coef.var 
##    1.4893271    0.0522833 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.3328 0.5739
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.0219 0.021895   12.71 0.00346 **
## Residuals   13 0.0224 0.001723                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1  20.21  20.212   5.803 0.0315 *
## Residuals   13  45.28   3.483                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value  Pr(>F)   
## df[[1]]      1 8.345e-05 8.345e-05   15.39 0.00175 **
## Residuals   13 7.051e-05 5.420e-06                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.021895 0.021895  12.096 0.00456 **
## df[[3]]      1 0.000679 0.000679   0.375 0.55179   
## Residuals   12 0.021721 0.001810                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.038963  1 21.5255 0.0005706 ***
## df[[1]]     0.011786  1  6.5114 0.0253839 *  
## df[[3]]     0.000679  1  0.3749 0.5517931    
## Residuals   0.021721 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.0253839025719477"
# cohort2.2_RQ_M_9 is probably an outlier
# dark_pt2 %<>% select(-cohort2.2_RQ_M_9)
# covar_SPF_pt1 <- c(cohort2_1_covar[[1]][c(6,8,9)], cohort3_1_covar[[1]][c(4,5,6)])

cycle = "Dark"
dark_pt2.RQ.rcf <- rcf_pipe(dark_pt2, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.6168776 1.0818660
## 
## $cut
## [1] 0.009299768

## [1] "SPF EC50: 0.902146907095478"
## [1] "GF EC50: 0.832689248018978"
t <- dark_pt2.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt1,covar_GF_pt1)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.2_RQ_M_9 = 0.67752, p-value = 0.01046
## alternative hypothesis: lowest value 0.84275526491303 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.2_RQ_M_9 = 0.073739, p-value = 0.3761
## alternative hypothesis: lowest value 0.744699285944684 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.74009, p-value = 0.009986
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.85735, p-value = 0.113
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 0.2754, num df = 6, denom df = 7, p-value = 0.1371
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.05380359 1.56852627
## sample estimates:
## ratio of variances 
##          0.2753989
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 3.6963, df = 13, p-value = 0.002689
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.0313418 0.1195076
## sample estimates:
## mean of x mean of y 
## 0.8969731 0.8215484
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7522 -1.2915 -0.6480  0.9957  4.7798 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   42.167      8.809   4.787 0.000355 ***
## df[[2]]      -14.521     10.263  -1.415 0.180605    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.089 on 13 degrees of freedom
## Multiple R-squared:  0.1334, Adjusted R-squared:  0.06679 
## F-statistic: 2.002 on 1 and 13 DF,  p-value: 0.1806
## 
## y = 42 + -15x, r^2 = 0.133

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  8.000000000  0.000000000  0.000000000  0.744699286  0.871081793  0.126382507 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  6.572387075  0.834724603  0.821548384  0.017086348  0.040402793  0.002335546 
##      std.dev     coef.var 
##  0.048327490  0.058824886 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 7.0000000000 0.0000000000 0.0000000000 0.8427552649 0.9241419801 0.0813867152 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 6.2788116371 0.9031934811 0.8969730910 0.0095857548 0.0234554971 0.0006432069 
##      std.dev     coef.var 
## 0.0253615234 0.0282745644 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   8.00000000   0.00000000   0.00000000  29.00000000  35.00000000   6.00000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 246.50000000  30.00000000  30.81250000   0.75555975   1.78661491   4.56696429 
##      std.dev     coef.var 
##   2.13704569   0.06935645 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##    7.0000000    0.0000000    0.0000000   26.3000000   30.4000000    4.1000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  199.4000000   28.1000000   28.4857143    0.5629127    1.3773978    2.2180952 
##      std.dev     coef.var 
##    1.4893271    0.0522833 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   2.815 0.1172
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df  Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.02124 0.021239   13.66 0.00269 **
## Residuals   13 0.02021 0.001554                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1  20.21  20.212   5.803 0.0315 *
## Residuals   13  45.28   3.483                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## df[[1]]      1 8.558e-05 8.558e-05   18.97 0.000779 ***
## Residuals   13 5.865e-05 4.510e-06                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq  Mean Sq F value  Pr(>F)   
## df[[1]]      1 0.021239 0.021239  12.651 0.00395 **
## df[[3]]      1 0.000063 0.000063   0.037 0.84979   
## Residuals   12 0.020145 0.001679                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##                Sum Sq Df F value   Pr(>F)   
## (Intercept) 0.0292315  1 17.4125 0.001293 **
## df[[1]]     0.0157706  1  9.3941 0.009808 **
## df[[3]]     0.0000629  1  0.0374 0.849794   
## Residuals   0.0201452 12                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
## [1] "Varies significantly with DV, p = 0.00980775246066405"

1 month after HF diet

# pt3 RQ ####

# cohort2.3_RQ_M_1 is an outlier.
total_pt3 %<>% select(-cohort2.3_RQ_M_1)
light_pt3 %<>% select(-cohort2.3_RQ_M_1)
dark_pt3 %<>% select(-cohort2.3_RQ_M_1)
covar_SPF_pt3 <- c(cohort1_3_covar[[1]][c(2,3)], cohort2_3_covar[[1]][c(4,5,6,7)], cohort3_3_covar[[1]][c(4,5,6)])
covar_GF_pt3 <- c(cohort1_3_covar[[1]][1], cohort2_3_covar[[1]][c(2,3)], cohort3_3_covar[[1]][c(1,2,3)])

pt = "1 month after HF diet"
cycle = "Total"
total_pt3.RQ.rcf <- rcf_pipe(total_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5034941 1.3791840
## 
## $cut
## [1] 0.0175138

## [1] "SPF EC50: 0.843815968423415"
## [1] "GF EC50: 0.814793610952982"
t <- total_pt3.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.3_RQ_M_12 = 0.3487, p-value = 0.4091
## alternative hypothesis: lowest value 0.785868886334165 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.3_RQ_M_2 = 0.55225, p-value = 0.1078
## alternative hypothesis: lowest value 0.779598208838337 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.9508, p-value = 0.6989
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.90215, p-value = 0.3868
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 2.1499, num df = 8, denom df = 5, p-value = 0.4148
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.3181729 10.3568953
## sample estimates:
## ratio of variances 
##           2.149949
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 1.8393, df = 13, p-value = 0.08881
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.004607271  0.057401018
## sample estimates:
## mean of x mean of y 
## 0.8425080 0.8161111
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.571 -2.741 -1.206  2.654  5.123 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    14.47      23.36   0.620    0.546
## df[[2]]        21.71      28.07   0.773    0.453
## 
## Residual standard error: 3.094 on 13 degrees of freedom
## Multiple R-squared:  0.04398,    Adjusted R-squared:  -0.02956 
## F-statistic: 0.5981 on 1 and 13 DF,  p-value: 0.4531
## 
## y = 14 + 22x, r^2 = 0.044

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 6.0000000000 0.0000000000 0.0000000000 0.7795982088 0.8372413178 0.0576431090 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 4.8966667261 0.8173628859 0.8161111210 0.0085068096 0.0218674502 0.0004341949 
##      std.dev     coef.var 
## 0.0208373428 0.0255324824 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 9.0000000000 0.0000000000 0.0000000000 0.7858688863 0.8876941418 0.1018252554 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 7.5825719507 0.8499652336 0.8425079945 0.0101843925 0.0234852513 0.0009334967 
##      std.dev     coef.var 
## 0.0305531776 0.0362645551 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   6.00000000   0.00000000   0.00000000  28.70000000  33.20000000   4.50000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 179.60000000  29.35000000  29.93333333   0.67806915   1.74303225   2.75866667 
##      std.dev     coef.var 
##   1.66092344   0.05548742 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  30.90000000  38.30000000   7.40000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 308.40000000  34.90000000  34.26666667   0.82276634   1.89730257   6.09250000 
##      std.dev     coef.var 
##   2.46829901   0.07203207 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.5224 0.4826
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.002508 0.0025085   3.383 0.0888 .
## Residuals   13 0.009639 0.0007415                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  67.60   67.60   14.05 0.00243 **
## Residuals   13  62.53    4.81                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 2.523e-05 2.523e-05   6.242 0.0267 *
## Residuals   13 5.254e-05 4.042e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.002508 0.0025085   3.241  0.097 .
## df[[3]]      1 0.000351 0.0003508   0.453  0.514  
## Residuals   12 0.009288 0.0007740                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.054279  1 70.1275 2.347e-06 ***
## df[[1]]     0.002325  1  3.0038    0.1087    
## df[[3]]     0.000351  1  0.4532    0.5136    
## Residuals   0.009288 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)


cycle = "Light"
light_pt3.RQ.rcf <- rcf_pipe(light_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5164798 1.3690690
## 
## $cut
## [1] 0.01705178

## [1] "SPF EC50: 0.796725491719892"
## [1] "GF EC50: 0.750865595089389"
t <- light_pt3.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.3_RQ_M_12 = 0.21909, p-value = 0.8691
## alternative hypothesis: lowest value 0.733507805063572 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort2.3_RQ_M_3 = 0.11486, p-value = 0.5711
## alternative hypothesis: highest value 0.821379469025812 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.98937, p-value = 0.9952
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.90861, p-value = 0.4273
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 0.76136, num df = 8, denom df = 5, p-value = 0.6956
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1126738 3.6676638
## sample estimates:
## ratio of variances 
##          0.7613565
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 1.666, df = 13, p-value = 0.1196
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01087743  0.08418140
## sample estimates:
## mean of x mean of y 
## 0.7969794 0.7603274
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.423 -2.636 -1.206  2.447  5.490 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    20.72      14.59   1.421    0.179
## df[[2]]        15.10      18.62   0.811    0.432
## 
## Residual standard error: 3.087 on 13 degrees of freedom
## Multiple R-squared:  0.04815,    Adjusted R-squared:  -0.02507 
## F-statistic: 0.6576 on 1 and 13 DF,  p-value: 0.432
## 
## y = 21 + 15x, r^2 = 0.0481

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  6.000000000  0.000000000  0.000000000  0.702764067  0.821379469  0.118615402 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  4.561964193  0.746089137  0.760327366  0.018450130  0.047427568  0.002042444 
##      std.dev     coef.var 
##  0.045193403  0.059439401 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  9.000000000  0.000000000  0.000000000  0.733507805  0.860417500  0.126909695 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  7.172814154  0.800611733  0.796979350  0.013144613  0.030311532  0.001555028 
##      std.dev     coef.var 
##  0.039433839  0.049479123 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   6.00000000   0.00000000   0.00000000  28.70000000  33.20000000   4.50000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 179.60000000  29.35000000  29.93333333   0.67806915   1.74303225   2.75866667 
##      std.dev     coef.var 
##   1.66092344   0.05548742 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  30.90000000  38.30000000   7.40000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 308.40000000  34.90000000  34.26666667   0.82276634   1.89730257   6.09250000 
##      std.dev     coef.var 
##   2.46829901   0.07203207 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.0367 0.8511
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df   Sum Sq  Mean Sq F value Pr(>F)
## df[[1]]      1 0.004836 0.004836   2.775   0.12
## Residuals   13 0.022652 0.001742               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  67.60   67.60   14.05 0.00243 **
## Residuals   13  62.53    4.81                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 1.648e-05 1.648e-05   3.435 0.0866 .
## Residuals   13 6.236e-05 4.797e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq  Mean Sq F value Pr(>F)
## df[[1]]      1 0.004836 0.004836   2.607  0.132
## df[[3]]      1 0.000393 0.000393   0.212  0.654
## Residuals   12 0.022259 0.001855               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.048143  1 25.9538 0.0002642 ***
## df[[1]]     0.003906  1  2.1055 0.1724047    
## df[[3]]     0.000393  1  0.2119 0.6535429    
## Residuals   0.022259 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)

cycle = "Dark"
dark_pt3.RQ.rcf <- rcf_pipe(dark_pt3, cycle, channel)
## [1] "Range:"
## $range
## [1] 0.5034941 1.3791840
## 
## $cut
## [1] 0.0175138

## [1] "SPF EC50: 0.871701952811977"
## [1] "GF EC50: 0.866975323021419"
t <- dark_pt3.RQ.rcf %>% lapply(EC50)
df <- data.frame(genotype=c(rep("SPF", length(t[[1]])), rep("GF", length(t[[2]]))), rcf=c(t[[1]], t[[2]]), covar=c(covar_SPF_pt3,covar_GF_pt3)) %>% set_rownames(NULL)

ggplot(df) + aes(x="", y=rcf, fill=genotype) + geom_boxplot() + ggtitle(paste(cycle, channel, "RCF EC50",  pt, sep = " "))

dixon.test(t[[1]])
## 
##  Dixon test for outliers
## 
## data:  t[[1]]
## Q.cohort2.3_RQ_M_10 = 0.041949, p-value = 0.25
## alternative hypothesis: lowest value 0.831931716444821 is an outlier
dixon.test(t[[2]])
## 
##  Dixon test for outliers
## 
## data:  t[[2]]
## Q.cohort3.3_RQ_M_1 = 0.44814, p-value = 0.2602
## alternative hypothesis: lowest value 0.805485833437532 is an outlier
f <- ggarrange(
  ggdensity(df[,1:2], x="rcf", color="genotype", palette=colors, add="mean", rug=TRUE),
  ggqqplot(df[,1:2], x="rcf", color="genotype", palette=colors)
  
)
annotate_figure(f, top = text_grob(paste(cycle, channel, "rcf distribution", pt, sep = " "), face = "bold"))

shapiro.test(t[[1]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[1]]
## W = 0.90685, p-value = 0.2944
shapiro.test(t[[2]])
## 
##  Shapiro-Wilk normality test
## 
## data:  t[[2]]
## W = 0.92123, p-value = 0.5142
var.test(t[[1]], t[[2]])
## 
##  F test to compare two variances
## 
## data:  t[[1]] and t[[2]]
## F = 0.53355, num df = 8, denom df = 5, p-value = 0.4096
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.07896096 2.57027073
## sample estimates:
## ratio of variances 
##          0.5335528
t.test(t[[1]], t[[2]], var.equal = T)
## 
##  Two Sample t-test
## 
## data:  t[[1]] and t[[2]]
## t = 0.40859, df = 13, p-value = 0.6895
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02685763  0.03938625
## sample estimates:
## mean of x mean of y 
## 0.8704050 0.8641407
res <- ANCOVA(df)
## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9276 -2.7650 -0.7837  2.5419  5.8042 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    33.84      26.02    1.30    0.216
## df[[2]]        -1.51      29.97   -0.05    0.961
## 
## Residual standard error: 3.164 on 13 degrees of freedom
## Multiple R-squared:  0.0001953,  Adjusted R-squared:  -0.07671 
## F-statistic: 0.002539 on 1 and 13 DF,  p-value: 0.9606
## 
## y = 34 + -1.5x, r^2 = 0.000195

## 
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## Some stats on the variables
## DV - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  6.000000000  0.000000000  0.000000000  0.805485833  0.899655802  0.094169968 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  5.184844213  0.871920663  0.864140702  0.014064760  0.036154616  0.001186905 
##      std.dev     coef.var 
##  0.034451485  0.039867911 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 9.0000000000 0.0000000000 0.0000000000 0.8319317164 0.9007889218 0.0688572053 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 7.8336450909 0.8710080145 0.8704050101 0.0083883275 0.0193435180 0.0006332764 
##      std.dev     coef.var 
## 0.0251649826 0.0289118081 
## Covar - IV:
## df[[1]]: GF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   6.00000000   0.00000000   0.00000000  28.70000000  33.20000000   4.50000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 179.60000000  29.35000000  29.93333333   0.67806915   1.74303225   2.75866667 
##      std.dev     coef.var 
##   1.66092344   0.05548742 
## ------------------------------------------------------------ 
## df[[1]]: SPF
##      nbr.val     nbr.null       nbr.na          min          max        range 
##   9.00000000   0.00000000   0.00000000  30.90000000  38.30000000   7.40000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 308.40000000  34.90000000  34.26666667   0.82276634   1.89730257   6.09250000 
##      std.dev     coef.var 
##   2.46829901   0.07203207 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.5065 0.4892
##       13               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000141 0.0001413   0.167  0.689
## Residuals   13 0.011001 0.0008462               
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## df[[1]]      1  67.60   67.60   14.05 0.00243 **
## Residuals   13  62.53    4.81                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value  Pr(>F)   
## df[[1]]      1 4.198e-05 4.198e-05   9.805 0.00795 **
## Residuals   13 5.566e-05 4.280e-06                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## <U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015><U+2015>
## ANCOVA
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000141 0.0001413   0.157  0.699
## df[[3]]      1 0.000210 0.0002098   0.233  0.638
## Residuals   12 0.010791 0.0008992               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.058262  1 64.7902 3.531e-06 ***
## df[[1]]     0.000349  1  0.3880    0.5450    
## df[[3]]     0.000210  1  0.2334    0.6377    
## Residuals   0.010791 12                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_ANCOVA(res)
 

By Sumeed Yoyo Manzoor

smanzoor@uchicago.edu



Back to top